using Plots

Efield = 5e8 # V/m
d=0.34e-9
thlist=[1,1.1,1.57,1.76, 1.95, 1.99]*pi/180

gr(reuse=true, size=(1000,1000))
h=plot(x->(-Efield*d*x))

fp=open("ScreeningCalc.txt", "w")

for Nlayer=2:6
    th=thlist[Nlayer]
    bw = 20e-3 # meV
    density = 8/sqrt(3)*(th/0.246e-9)^2
    @show Cq = density/bw*(1.602e-19)/Nlayer

    Cg = 8.85e-12 / d

    g = Cg/Cq # Quantum capacitance ratio
    @show V0 = 8.85e-12 * Efield / Cg
    B = zeros(Nlayer,1); B[1] = V0*g; B[end] = -V0*g
    A = diagm(0 => [1+g;ones(Nlayer-2)*(1+2g);1+g], 1=>ones(Nlayer-1)*(-g), -1=>ones(Nlayer-1)*(-g))
    V = (A\B)[:]

    z=-(Nlayer-1)/2 : (Nlayer-1)/2
    plot!(z, V,marker=(:circle,5), ylims=(-0.1,0.1))
    for i=1:Nlayer
        println(fp, "$(z[i]) $(V[i]) $(Nlayer)")
    end
    println(fp)
end
close(fp)
h
